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Abstract We have formulated and experimentally demonstrated an improved algorithm for design of 
arbitrary two-dimensional holographic traps for ultracold atoms. Our method builds on the best previously 
available algorithm, MRAF, and improves on it in two ways. First, it allows for creation of holographic 
atom traps with a well defined background potential. Second, we experimentally show that for creating 
trapping potentials free of fringing artifacts it is important to go beyond the Fourier approximation in 
modelling light propagation. To this end, we incorporate full Helmholtz propagation into our calculations. 



Optical dipole traps have become ubiquitous in disciplines ranging from live cell manipulation in biophysics [T| to 
single atom manipulation for quantum information processing |2j. In the field of ultracold atomic gases a variety of 
optical potential shapes, of ever increasing complexity, are used for fundamental studies of many-body physics [3] in 
different optical crystal lattices [H [5] , reduced dimensionality |6', , and non-trivial trapping topologies [H |9] . Optical 
sculpting is also likely to offer great benefits for trapped-atom interferometry with Bose-Einstein condensates (BECs) 
[TUj and for creating "atomtronic" [TT] optical circuits. 

The simplest optical trap for an ultracold atomic gas is formed by a single focused Gaussian laser beam [T^ . Natural 
extensions include standing- wave optical lattice potentials produced by the interference of multiple laser beams [4], 
and ring traps produced using Laguerre- Gauss laser modes |51 [H]- However there is so far no universal approach 
for creating an arbitrary optical trapping potential on the required micrometer scale. Existing techniques fall into 
three categories: (1) a time-averaged potential can be created by the fast scanning of a focused laser beam [T3], (2) 
an intensity mask or micromirror-device can be imaged onto the trapping plane of the atoms |141 115j . and (3) in 
the holographic method one manipulates the phase of the laser beam so as to create the desired intensity pattern 
after further propagation of the light. None of these methods is currently a clear overall winner; their relative merits 
depend in practice on various factors, such as the desired spatial resolution, efficiency of use of laser light, and temporal 
stability of the light pattern. 

In this paper, we provide a computational and experimental procedure for generating arbitrary potentials based on 
an improved holographic, phase only method. We illustrate our technique by producing high fidelity patterns which 
demonstrate suitable characteristics for atom trapping, opening a versatile toolbox for optical manipulation in exotic 
geometries. 

As shown in Figjlja), we use a spatial light modulator (SLM) to imprint a custom phase pattern (kinoform) on a 
laser beam in order to form a diffraction pattern of the desired shape in the trapping plane. Since this is a diffractive 
method, the kinoform steers the light into required trap shapes without blocking out any of the incident power, making 




FIG. 1. (a) Illustration of the optical system used for creating the traps. CoUimated laser light is reflected from the SLM 
through a Fourier 2/ arrangement to form the intensity pattern in the trapping plane. / denotes the lens' focal length, see the 
Methods section for definition of the remaining symbols. The pattern shown in the trapping plane is the atomtronic OR-gate, 
which we use to test our method, (b) A cartoon of the constraints and the variables to be solved for. The modulus of the light 
field is fixed in the SLM plane, and we wish to solve for the kinoform to produce the correct modulus with no constraints on 
the phase in the trapping plane. 
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the method potentially very efficient. In some cases, the required kinoform can be calculated analytically [Slin], but in 
general, this is a highly non-trivial numerical problem. The difficulty arises from the fact that we have incomplete a 
priori knowledge of the light field in both planes: In the trapping plane only the modulus of the light field is relevant 
for optical trapping, so we place no constraint on the phase; similarly, at the SLM, we impose the unaltered Gaussian 
modulus of the incident laser profile, leaving the phase as the variable to be calculated [Fig. [lib)]. 




The numerical task of calculating the kinoform is approached as an optimisation procedure which can be achieved 
by "steepest descent" methods such as direct binary searches [T6U18) . genetic algorithms [19] or Gerchberg and 
Saxton's iterative Fourier transform methods [20] • Our method is based on the Gerchberg-Saxton algorithm (see 
the Methods section) since it is very computationally efficient and is easily executed on a graphics processing unit. 
So far, the best computational results were obtained by Pasienski and DeMarco [21], using an adaptation of the 
Gerchberg-Saxton scheme known as the Mixed-Region-Amplitude-Freedom (MRAF) algorithm. The idea of MRAF 
is to enhance convergence of the iterative algorithm in one region of the trapping plane by giving up control of the 
remaining regions |21H24j . This approach produced excellent computational results for a range of generic trap shapes, 
and some shapes have also been successfully realised with laser light [25] [26] . However, one practical limitation of 
MRAF is that in order to achieve good convergence in bright areas of the trapping plane (green areas in Fig. T]), we 
must surrender control over almost all of the dark, background regions. In particular, the original MRAF calculations 
provide for only a very narrow region defining the zero background potential to which all the potentials are referenced. 
This could lead to experimental problems such as inefficient loading of atoms into, or percolation of atoms out of 
the trap; this issue is particularly important if the desired trap has sharp edges so the region of uncontrolled light 
intensity is immediately adjacent to the region of high atomic density. If the region of zero potential is extended, the 
fidelity of MRAF kinoforms is severely compromised. The first result of this paper is to show that this problem can 
be eliminated by properly redefining the target potential so as to include both the conventionally defined trapping 
region and a sufficiently large well defined background region (a "canvas" on which an arbitrary potential landscape 
is "drawn"). We formulate an offset-MRAF (OMRAF) algorithm through a simple set of practical rules which obey 
the limits imposed by Nyquist's theorem and minimise the occurrence of convergence-stalling optical vortices. 

The second result of this paper is to demonstrate that the paraxial approximation often employed in literature 
[HJ [57] can lead to undesirable artifacts in the traps. We illustrate the improvement that can be made by direct 
numerical calculation of the light propagation using the Helmholtz equation. Using our computational algorithms 
and a simple apparatus, we experimentally demonstrate examples of optical traps suitable for atomic experiments, 
and compare the results to previous realisations [25, 26]. 

The paper is organised as follows: First we compare computational results produced using the existing MRAF and 
our OMRAF algorithms. Then we present and analyse experimental images of light fields shaped by our method. 
Throughout this investigation, we primarily test our technique by creating an atomtronic OR-gate trap shape jllLI21l 
[28] (see Fig[T]). To illustrate the universality of our method we also include the final results for a uniform square trap 
and an annular BEG stirrer [5TJ |^ . 



Within the Gerchberg-Saxton framework, the core algorithm relies on linking the moduli in the SLM and trapping 
planes by simulating the light propagation back and fourth between these planes. After each propagation, we manually 
impose the known modulus constraints while leaving the phase to converge on the required solution. The mathematical 
details are presented in the Methods section, and here we just briefiy review the algorithms for kinoform calculation 
before presenting the numerical and experimental results. 

It is well documented that convergence of the naive Gerchberg-Saxton algorithm is highly erratic because the 
modulus constraints in each plane are non-convex [29] [30]. The MRAF algorithm improves on this by defining 
a "drawing" V in the trapping plane containing all the points in which the desired potential is non-zero, and an 
additional "canvas" region, C, of zero potential around the drawing [Fig. [2ja)]. The modulus constraint is then 
manually imposed only inside C and V at each iterative step, thus reducing the number of constraints the algorithm 
aims to satisfy. 

In the original MRAF paper, region C is limited to a tight border around V which has a typical width of around 
1% of the trapping region length scales. Although C is a featureless region of zero potential, it is an integral part 
of the trap since it defines the zero to which all the pixels in V are referenced. For practical purposes such a small 
canvas can be a severe limitation. Our method aims to expand C to approach its theoretical limit which is set by 
the Nyquist theorem to be 25% of the area of the simulated trapping plane (see the Methods section) . In this large 
canvas regime, MRAF convergence stagnates and produces poor solutions even after many iterations (see Fig. [21 




RESULTS 



The Gerchberg-Saxton algorithm and MRAF 
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FIG. 2. Computational results for traps with a well defined zero potential background, (a) Reconstructions of the OR-gate 
trap with MRAF and OMRAF methods. We also show details of the local reconstruction error (bottom insets) and the local 
phase (right insets). In the phase plots, obtained after 30 iterations, we use black to white to represent to 2t7 and red dots to 
mark vortex cores, (b) Illustrates the method for picking the optimal offset. A, measured in units of the peak amplitude in the 
target pattern. We perform a single iteration of the algorithm, for different offsets, to produce patterns such as those shown in 
the insets. We then choose the lowest value of A (indicated by the dashed line) which gives no vortices within C or D in this 
first iteration, (c) Improvement in the trap reconstruction associated with the offset. 



Optical vortices and OMRAF 

We observe the problems of the MRAF algorithm for extended C to be due to the formation of large populations 
of optical vortices during the early iterations (approximately 1 vortex for every 10 pixels in Fig. [2|a)). Optical 
vortices are topological features of the light field corresponding to a phase winding around a point of zero intensity. 
Large vortex populations are problematic because individual vortices are difficult to eliminate in further iterations of 
the algorithm. The reason for this is that the Gerchberg-Saxton algorithm monotonically reduces the RMS error in 
each iteration |31j . but topological (un) winding operations cause disruption across the entire kinoform. Such global 
disruption is not in keeping with the monotonic error reduction in later iterations |32j . Therefore any erroneous 
vortices are "frozen in" early on, and the algorithm gets stuck in a local RMS minimum corresponding to a particular 
vortex distribution. 

One possibility for handling vortices is by pairwise creation and annihilation, which requires only local operations; 
computational algorithms which encourage such pairwise operations are given in [33. However, we consider a much 
simpler solution. Our method simply offsets the trapping pattern inside C and I? by a uniform intensity, |Ap, to 
remove all points of zero light intensity. This redefines the zero of the potential but does not change the physics of the 
trap. With an offset intensity, we know a priori that we should aim to create no vorticies at all inside C, so we feed 
the algorithm with an initial kinoform which contains no phase winding. Specifically, we choose a parabolic initial 
kinoform [21, 24 which defocuscs the light into a patch with a characteristic size set by C. 

The optimal value of the offset intensity is calculated semi-empirically as follows: After just one iteration, the 
intensity in the trapping plane adopts approximately the correct shape, but exhibits considerable fluctuations around 
the desired intensity [Fig. |2jb)]. If any of these fluctuations bring the intensity locally to zero, a vortex may form at 
this point. Therefore, we increase |Ap until the trapping plane contains no vortices inside C after the first iteration 
[Fig. [2jb) and (c)], and then trust that very few vortices will be formed in subsequent iterations. For maximum 
light-usage efficiency, we want to minimise the proportion of the incident light which is steered into the featureless 
background, i.e. choose the minimal |A| which gives satisfactory results. We find that for all the patterns we tested a 
suitable offset was | A| '-^ 10 — 15% of the maximum amplitude in the trapping plane (i.e. | A|^ 1% of the maximum 
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intensity). The computational simulation of the trapping potential produced using the OMRAF algorithm is shown 
in Fig. [2j Our method reduces the number of vortices seen in C after 30 iterations by a factor of ~ 20 compared to 
MRAF. This allows excellent convergence of the algorithm, leaving only 4% RMS error after 30 iterations and 1% 
RMS after 1000. 

One compromise associated with the OMRAF method is a reduction in the efficiency of light use. We define 
efficiency as the ratio of the integral of the trap intensity above the background offset to the total integral across the 
trapping plane. Including the non-zero A causes a drop in efficiency by a factor of approximately 2, from 43% to 
24%, for the OR-gate. Nevertheless, our method remains more efficient than intensity masking methods such as |14j . 
which report 3% efficiency for simple patterns. We calculate that for more complex patterns such as the OR-gate, 
intensity masking would be less than 1% efficient. 



Experimental Realisation witli Laser Liglit 

We test our OMRAF algorithm experimentally by producing light patterns using 532 nm laser light in the arrange- 
ment shown in Fig. [ija). We use the standard SLM-based Shack-Hartmann algorithm [53] to crudely correct for 
the low spatial frequency phase aberrations in the optical system and characterise our input laser beam. We then 
use an active feedback algorithm j25j to optimise the final experimental patterns. This feedback routine is essential 
to remove the "sine envelope" caused by the pixellation of the SLM, which would otherwise globally modulate the 
intensity pattern. Finally, in practice, we find that when we try to produce experimental patterns with C covering the 
maximum theoretically permitted area (25% of the trapping plane), we cannot achieve patterns with less than 20% 
RMS fluctuations. This is significantly improved in the results presented below by reducing C to cover only 10% of 
the trapping plane. 

In the previous section, we presented computational results generated under the paraxial approximation in which 
the propagator mapping the light in the SLM plane to the trapping plane is given by a scaled Fourier transform 27] . 
Using this approach in our experiments with an f — 200 mm lens, we are able to produce a trapping pattern with RMS 
error of « 11% [Fig. [3]ja)]. As highlighted in Fig. |3][b), the most prominent form of error is semi-regular fringing. We 
suggest that the source of this error is the inadequacy of the paraxial approximation, and use a numerical Helmholtz 
solver (see the Methods section) to test this hypothesis. Speciflcally, we use the OMRAF algorithm under the 
paraxial approximation to produce a "Fourier-generated" kinoform, and then input this kinoform into our numerical 
Helmholtz propagator to simulate the experimental apparatus. This indeed produces the fringing pattern similar to 
the one observed experimentally, which gives us confidence to continue with the Helmholtz method. We thus replace 
all instances of the Fourier propagator in the OMRAF algorithm with a Helmholtz propagator. 

In this way, we are able to eliminate the erroneous fringing, and reduce the RMS error of the experimental pattern 
to 7% [see Fig. ^a)]. As shown in Fig. [sfc), we can produce other trap shapes with similar fidelity. Thus, for 
traps with sharp edges and an extended canvas region, we achieve RMS variation comparable to the 4% fluctuations 
previously seen only for simple smooth traps ;25j. In Ref. |3S], it has already been shown that 5% RMS errors are 
sufficiently low for atomtronic applications. We therefore believe that our methods show great promise for future 
applications. 



DISCUSSION 

We have presented a simple method for producing optical traps of arbitrary 2D profile suitable for use in cold atom 
experiments. A useful optical potential must be referenced to a well defined background, and we showed that in order 
to do this, we must modify the existing algorithms to include a non-zero background light intensity. Specifically, 
we developed the OMRAF algorithm guided by the principle that points of zero intensity seed convergence-stalling 
vortices, and therefore should be avoided. In computer simulations we produce traps covering 25% of the trapping 
plane (the Nyquist limit) with 4% RMS error after only 30 iterations. We also took the next step of realising these 2D 
profiles in laser light. Here, we remove systematic aberrations using a Shack-Hartmann method and active feedback. 
In addition, we found that using a beyond-paraxial Helmholtz solver improves the fidelity of the experimentally 
produced light patterns. 
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FIG. 3. Experimental results with 532nm laser light, (a) Experimental patterns based on four different computational methods 
(for display we crop the trapping plane around region C). We observe low quality traps with the MRAF algorithm, and we are 
unable to improve them by adjusting the propagator (indeed, in this high RMS error regime, the additional complexity of the 
Helmholtz propagator can be detrimental). However, incorporating a Helmholtz propagator into our OMRAF algorithm allows 
us to produce the OR-gate with only 7% RMS error, (b) Illustrates the need for a beyond-paraxial propagator. The left panel 
shows an experimental pattern based on OMRAF calculation under the paraxial approximation. The right panel shows that 
similar fringing is reproduced when the same paraxial kinoform is fed into a Helmholtz propagator to simulate the experiment, 
(c) To test the universality of our method, we also show experimental OMRAF patterns cropped to the boundary of C for a 
uniform square trap and an annular BEG stirrer. 



METHODS 



Modified Gerchberg-Saxton Algorithms 



Both the MRAF and the OMRAF approach have the Gerchberg-Saxton algorithm at their core. This algorithm 
attempts to obtain the kinoform, 0(r), which must be applied to the input Gaussian field, Eo{r), in the SLM plane 
(spanned by r) to obtain the desired trap shape, r(R), in the trapping plane (spanned by R) at z = 2/ [sec Fig. 
pTa)]. These two planes are related by a projection operator, so we may write the implicit equation for (j){r) as: 



|^[|£;o(r)|exp{*</)(r)}]| =r(R). (1) 

The Gerchberg-Saxton algorithm solves this equation iteratively as illustrated in Fig. |4j 

MRAF relaxes the modulus constraint outside the trapping region by defining a drawing, V (defined by r(R) > 0) 
together with a narrow canvas, C, of zero intensity around V, and then applying the following algorithm: 



^in) ^ f mT(R) exp {i arg (a^^^ (R)) } R e 2? U C 
' \ (l-m)4"^(R) R^ 2?UC ' 



We find that a mixing parameter m = 0.4 gives the minimum RMS deviation for all the patterns trialled (see also 

OMRAF further modifies this algorithm in the following ways: First we expand the canvas so that the trapping 
region covers the maximum theoretical area set by Nyquist's theorem (25% of the whole trapping plane - see the next 
section). Secondly, throughout C and V, we offset the desired light intensity according to: 



r(R) ^ ^/TCR)^ + A2 



(3) 
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FIG. 4. Diagrammatic representation of the core Gerchberg-Saxton algorithm. Refer to the text for definition of symbols. 



This shifts both the trapping potential and its zero reference level equally, and does not change the physics of the 
trap. 



Maximum Canvas Size 



In order to best use the degrees of freedom offered by an x pixel SLM, we should allow region C to cover as 
many pixels as possible in the discretized trapping plane. This theoretical limit is found by the following Nyquist 
argument: The N x N SLM plane can be fast-Fourier-transformed to an iV x iV trapping plane with a diffraction 
limited spot size of exactly 1 pixel. The Nyquist limit states that the diffraction limit should be reduced to half a 
pixel (the smallest feature in the trap shape). This can be done by padding the computer-simulated SLM plane with 
zeros to increase its size to 2N x 2N. This means that only 25% of the simulated SLM plane is actually covered by 
the physical SLM, so we can at best control only 25% of the trapping plane. 



A Computationally Efficient Helmholtz Solver 



We see improvement in the experimentally realised trapping intensities when the following method is used to model 
the propagation operator ^ in the OMRAF iterations: 

The Helmholtz equation for light propagation in free space gives that the Fourier transform of a light field, E{r; z), 
is modified by a phase factor Q{^; z) as it propagates in the +z direction [27] 

^ [E{v; z); ^] = ^ [£;(r; 0); exp {tkz^l-iM)^}, (4) 

^ ^ ^ 

e(C;2) 

where [g{r);^] denotes the Fourier transform (FT) of g, making the Fourier variable pairing r ^ explicit, and 
k = 2t:/X denotes the wave vector of the light. The routine for propagation of Eq via the apparatus in Fig. [ija) can 
then be conceptually decomposed into 3 stages: 

"Propagate FT from SLM to lens" -01 (0 
"Apply Lens in Fourier space" "02 (j?) 
"Propagate FT to trapping plane" ipsif]) 



^[i?o(r);e]e(C;/) 



tpi ® L 



MvMvJ) 



Then: ^Helmholtz [Eo{r);R] = [03(?y); R] , (5) 
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where ® represents convolution, and L represents the FT of the phase pattern imprinted by the lens, which, after 
aberration correction, we assume to be: 



(6) 
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